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Abstract. The diffusion-limited flow of hydrogen in the thermosphere is obtained by a I ~~ 
Monte Carlo calculation. The resulting density distribution has a steeper gradient than cur- f £ If- 
rent theories predict, even taking omission of gravity into account. The departure is ascribed - 

to increasing nonvalidity of the Chapman-Enskog diffusion coefficient as density changes over 
a mean free path become significant. 


Introduction 

The altitude distribution of minor nonperma- 
nent constituents of the atmosphere and their 
rate of escape have been treated by several 
authors, among them Spitzer [1952], Nicolet 
[1960], and Mange [1961]. In general these 
authors base their analysis on the hypothesis 
that a minor constituent exists in hydrostatic 
equilibrium under its own partial pressure. 
Mange [1961] and Bates and Patterson [1961] 
have taken flow into account to arrive at a 
steeper density gradient for hydrogen than 
would exist without flow, but hydrostatic equi- 
librium is still assumed. In an earlier study, not 
based on this hypothesis, Helg e-Petersen [1928] 
concluded that flowing helium did not attain 
hydrostatic equilibrium and therefore would 
decrease in density more rapidly with height 
than if such equilibrium were attained (see also 
Mitra [1952, p. 24]). 

Recently some indirect evidence has been 
obtained for the hypothesis that the hydrogen 
density gradient may be steeper than the 
current literature predicts. The intensity of 
night-sky Lymans radiation was found in a 
rocket flight (J. M. Coogan, private communi- 
cation) to decrease quite rapidly above 200 
km, by 30 per cent in 50 km. Such a decrease 
may reflect a steep density gradient, though the 
radiative transfer process is a complicated one. 

Donahue and Thomas [1963], also, suggest 
that the Lyman-** radiation appears to origi- 
nate, in part, by single scattering from an ex- 
tensive hydrogen envelope which is denser than 
would be expected on the basis of current esti- 
mates of the escaping flux. 

These factors motivated the present work, 


in which the Monte Carlo method is applied to 
the problem of hydrogen diffusion in the ther- 
mosphere, the object being to avoid hypotheses 
such as that of hydrostatic equilibrium and to 
simulate actual collisions as closely as possible. 

The physical model dealt with here is prob- 
ably best described by a comparison with the 
analytic treatments of Nicolet [1960] and 
others. In such treatments the vertical trans- 
port (flow) velocity of a minor constituent of 
the atmosphere is given by 


+ ( i + a T ) (1) 

where 

D = diffusion coefficient. 

ni = number density of minor constituent. 

mi — mass of minor constituent. 

ct T = thermal diffusion coefficient. 

The first term on the right-hand side repre- 
sents usual diffusion due to concentration gradi- 
ent. It is the major term and the one essentially 
calculated here. The second and third represent 
the effect of gravity, and the fourth, a relatively 
minor term, represents thermal diffusion. 

In the present calculation we do not include 
the effect of gravity; i.e., between collisions no 
force is assumed to act on the H atoms. There- 
fore the flow is limited only by diffusion. The 
effect of this admittedly severe simplification 
will be discussed below. 

The thermal diffusion term is believed to be 
automatically present in the calculation, since 
it results from exchange of momentum between 
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the H atoms and N 2 molecules during elastic 
collisions. 

Application of Monte Carlo Method 

We assume for simplicity that the major at- 
mospheric constituent is nitrogen, whose den- 
sity and temperature depend on altitude as in 
the Cospar atmosphere. The region dealt with 
lies between 100 and 500 km. This region is 
divided into 100 zones ; each zone is 4 km thick 
and has the appropriate Cospar density and 
temperature. 

Hydrogen atoms are introduced at levels of 
interest and made to undergo collisions with 
successive nitrogen molecules produced from 
appropriate Maxwell-Boltzmann (MB) distri- 
butions. In each collision, scattering is assumed 
to be isotropic in the center-of-mass system. 
The cross section for H-N z collisions is not 
available either experimentally or theoretically, 
so that an assumption about the interaction 
potential is necessary. A Lennard- Jones (6-12) 
potential was assumed here, as will be discussed 
below. Conveniently, the possibility of forma- 
tion of an H-N 2 complex may be ruled out be- 
cause of the high stability of N*. A detailed 
description of the process constructed follows. 

Each zone z has an MB distribution of N 3 
molecules determined by its temperature T. A 
hydrogen atom is injected into zone z — 1, 
with velocity chosen at random from the MB 
distribution of this zone. An N 2 molecule is pro- 
duced with velocity also chosen at random from 
the MB distribution. The collision cross section, 
which depends on the relative velocities, is com- 
puted, and it leads to an actual path traversed. 
The collision then occurs, and the resulting 
velocity of the H atom is obtained. A new N 2 
molecule is produced, and the process is re- 
peated. Transitions from one zone to the next 
are handled much as in Cashwell and Everett 
[1959]. However, in the upper zones, which 
have low density, a collision-free path may 
pass through many zones, so that a procedure 
was incorporated to add traversal times to the 
intervening zones. 

The Hydrogen-Nitrogen Collision Process 

(a) As was mentioned above, scattering is 
assumed to be isotropic in the center-of-mass 
system. The collision formulas of Cashwell and 
Everett [1959, pp. 56-60] were used in their 


general form, i.e. as applied to collisions with a 
moving target particle rather than a fixed 
target. 

( b ) A cross section for collisions between H 
atoms and N 2 molecules was obtained by means 
of 'combining rules/ which are discussed by 
Hirschf elder et at. [1954, p. 567]. The inter- 
molecular potential assumed was the Lennard- 
Jones 6-12 potential, 

4>(r) = 4«[(«r /r)’ 2 - (<r/r) 8 ] (2) 

where 

a = distance at which potential is zero. • 

€ = depth of potential minimum. 

For two dissimilar molecules, the parameters a 
and €, according to the combining rules, are to 
be taken as 

o’ 1 2 = $(01 + <7a) 

£12 = y / ' £1*2 

where <r lt e,, a 2) c 2 are the parameters that hold 
for like-molecular collisions of molecules 1 and 
2, respectively. 

The collision cross section for this intermolec- 
ular potential is (in the notation of Flirschfelder 
et at. [1954]) 

q (1) - We (l> * 

where the reduced cross section is a function 
only of the reduced kinetic energy of relative 
motion g * 2 : 

g* 2 = W s /« 12 

g being the initial relative speed of the colliding 
molecules and jjl the reduced mass. 

The quantity Q (1) * is plotted on page 558 of 



Fig. 1. Steady-state fluxes in two regions, 
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Hirschf elder et ah [1954]; values were taken 
directly off this curve for the purpose of the 
calculation. 

The parameters <s x = <s H , ei = e H , = cr N , } 
e 2 = e Na are thus needed. <r Ni and e Na are given 
in Table DA of Hirschf elder et ah [1954, p. 1111] 
as 3.75 X 10~ 8 cm and 1.10 X 10 -14 erg, respec- 
tively. For hydrogen the Lennard- Jones param- 
eters are not given, but the H% potential is well 
known and is given, for example, by Herzberg 
[1950]. On the assumption that this actual 
potential closely resembles a Lennard-Jones 
potential, the effective values of cr H and e a are 
0.39 X 10~ 8 cm and 7.16 X 10~ 12 erg, respec- 
tively. The second value is large because Hz is a 
tightly bound molecule, which means that 
ei 2 — (e A v e H ) in is correspondingly large, 
implying a very large cross section indeed for 
small relative velocities. This is considered 
reasonable, however, and in any case the above 
combining rules seem to be the only ones available 
for obtaining intermolecular potentials for such 
dissimilar molecules. 

(c) The effect of ignoring atomic oxygen 
may be briefly discussed at this point. The cross 
section for H-0 collisions must be large, since 
the OH radical can be formed in three-body 
collisions. Thus the presence of 0 will hinder 
diffusion of H, decreasing the density gradient 
somewhat. It would be feasible to actually in- 
clude a given proportion of 0 atoms in the 
Monte Carlo calculation to verify this quali- 
tative statement, but this was not done here. 
Bates and Nicolet [1950] estimate that the 
density of O may be as high as 10“ atoms/cm 3 
at 95 km, the O/N* ratio increasing with alti- 
tude, so that such a refinement might be worth 
while. 

Formation and dissociation of OH can prob- 
ably be neglected above 100 km, since, accord- 
ing to Bates and Nicolet [1950], its density is 
only 10 4 atoms/cm 8 at 90 km and drops rapidly. 

Calculation of Transmissions and Fluxes 

The Monte Carlo method can be used to ob- 
tain the density distribution of particles diffus- 
ing through an extended dense inhomogeneous 
medium, but it runs into the difficulty that in- 
sufficient numbers of particles penetrate the full 
extent of the medium. Kahn [1956] and others 
have devised methods of weighting (‘splitting/ 
etc.) to circumvent this difficulty. 



Fig. 2. Medium consisting of n regions. 

The alternative method used here describes 
the diffusion problem in terms of a transmission 
and reflection problem, utilizing the fact that 
transmission coefficients of finite regions of the 
medium are easily obtainable by the Monte 
Carlo method. The transmission relations have 
been essentially derived (in more complex form) 
in fields like multilayer optics [Heavens, 1955] . 

Consider two 'slab’ mediums in contact, whose 
transmission coefficients separately are 7\, T 2 ; 
i.e., out of a flux of Nv x particles cm" 2 sec" 1 
injected into the first slab, a flux T x Nv x will be 
transmitted. Of this, in turn, a fraction will 
be further transmitted through the second slab. 

There are multiply reflected and transmitted 
fluxes (Figure 1), and if the system is in a 
steady state (as much flux leaving the system 
as entering), we find, when these partial fluxes 
are totaled, that: (a) the transmitted flux at 
plane i = 3 is given by 

NviTiTz/il — RiR 2 ) 

where R x = 1 — T x ; R 2 = 1 — T 2 , so that the 
over-all transmission is given by T X TJ\ — R X R Z ; 
and (b) the forward flux at plane i — 2 is 

Nvi Ti/(l — RiR 2 ) 

The backward flux at plane % — 2 is 

NviT x R 2 /(l — R\ R 2 ) 

so that the flux difference at this plane is 


TABLE 1 


Region, km 

Transmission 

140-152 

0.0040 ± 0.0015 

152-168 

0.0074 d= 0.0026 

168-184 

0.0136 ± 0.0034 

184-200 

0.0183 dr 0.0037 

200-296 

0.0095 =b 0.0037 

296-496 

0.0231 ± 0.0061 




160 


S. O. KASTNER 



Fig. 3. Density distributions of assumed main nitrogen atmosphere (dashed line), Bates- 
Patterson hydrogen atmosphere (dash-dot line), hydrogen atmosphere derived here (solid 
line) ; k is a different (unspecified) parameter for each curve. 


where Ri = 1 — TV 

By continuing this iteration process, the over- i 
all transmission is obtained : 

T s T n _i 

The forward and backward fluxes at any 
given plane i can also be obtained by comput- 
ing the combined transmissions of the regions 
to the left and to the right of the plane and by 
applying the two-region formulas above. 

Let these fluxes at the plane i be denoted by 
7ii + Vi and nr v i} respectively. The total density 
at the plane i is (n» + + nr ) ; thus the total flux 
+ nr) must be divided by the mean 
particle velocity v it which is not obtainable from 
the transmission measurements alone. To find 
Vi we must relate the transmission description of 


NvxTxT */ (1 - B X B 2 ) 

in agreement with (a) . 

Now consider a slab medium of thickness l in 
the z direction. Divide it by planes into n re- 
gions, as in Figure 2, with the planes labeled by 
the index i. The transmission coefficients T lf 
T 2 , T n of the regions are assumed to be 
known (by injecting particles into one such 
region at a time). 

The combined transmission of the last two 
regions is 

Ti = T n -xTj(l - R n . 1 R n ) 

The combined transmission of the last three 
regions is then 

T 2 = T n _ 2 T 1 /(l - Bn- &x) 
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Altitude, 

km 

i 


Fr 

Si a (Ft* + Fr) 

Mi m (Ft* - Fr) 

140 

1 

1.0000 

0.9985 

1.9985 

0.0015 

152 

2 

0.6281 

0.6266 

1.255 

0.0015 

168 

3 

0.4278 

0.4263 

0.854 

0.0015 

184 

4 

0.3005 

0.2990 

0.600 

0.0015 

200 

5 

0.2204 

0.2189 

0.439 

0.0015 

296 

6 

0.06466 

0.06316 

0.128 

0.0015 

496 

7 

0.0015 

0.0 

0.0015 

0.0015 


the system to its description as a diffusing 
process. 

A given atom, moving in a random walk, 
undergoes a displacement £ along the z axis in a 
time At. £*, the mean square of such displace- 
ments for all atoms originating in a volume A 7, 
is related [. Kennard , 1938] to the diffusion coeffi- 
cient Z), by 

D t = lim £- (3) 

At- 0 AZ 

The mean displacement £,• at the plane t also 
gives the velocity v 4 required above; i.e., 


h = lim 

At-0 At 


( 4 ) 


The total density at plane t is then given by 
( n { + + n^)vj 


7li + + 7li = 


Vi 


MfVu Ta± *** > ^n) 


Vi 


(5) 


where /, is constructed as described above 
(combining transmissions to right and to left of 
the given plane, etc., to arrive at the sum of 
forward and backward fluxes). 

The mean time 1 4 — 1 /v 4 spent per atom in a 
unit volume at plane i is more directly gotten in 
practice than v if so that (5) is used in the form 


rii = I*(n< + + n 4 )vi (5a) 

In the actual computer program, atoms 
injected at each plane are always picked from an 
isotropic and Maxwellian distribution, a distribu- 
tion that cannot strictly hold in the actual 
diffusion process. The error resulting from this 
approximation should be small, however. 


Results 

(a) Transmissions and fluxes . Several thou- 
sand H atoms were successively injected at alti- 
tudes of 140, 152, 168, 184, 200, 296, 496 km, 
corresponding respectively to planes i = 1, 2, 3, 
4, 5, 6, 7. The resulting transmissions are given 
in Table 1. These results were used to find val- 
ues for forward and backward fluxes, which are 
given, with their sums and differences, in Table 2. 

The over-all transmission of the atmosphere 
from 140 to 500 km is therefore 0.0015, which 
is also the net upward flux when unit flux is 
injected at 140 km. (Extra significant figures 
included in the table entries are retained for 
calculational purposes and are not intended to 
represent unwarranted precision). 

( b ) Diffusion times and relative densities. 
At each of the planes i y H atoms were succes- 
sively injected and timed as they diffused 
through a 4-km distance. The mean times and 
resulting relative densities are presented in 
Table 3. 

It is to be noted that t 7 > U, which appears 
to be a real effect, explainable physically as due 
to freedom of the atoms to travel horizontally 
for long distances without collisions. 


TABLE 3 


Altitude, 

km 

i 

Time U To Diffuse 
4 kin, sec 

Relative 
Density USi 

140 

1 

68.1 

136 

152 

2 

25.5 

32 

168 

3 

16.4 

14 

184 

4 

11.5 

6.9 

200 

5 

9.0 

3.95 

296 

6 

2.1 

0.27 

496 

7 

4.4 

0.0066 
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i 

'dn(H) 

i) 

4 

D it 

cmVsec 

D, 

'dn(H) 

dz 

dz 

no. /cm 

1 

48.7 


1.54 X 10» 

7.5 X 10 10 

2 

17.8 


3.11 X 10» 

5.5 X 10“> 

3 

6.54 


5.35 X 10» 

3.5 X 10>° 

4 

2.8 


8.22 X 10» 

2.3 X 10 10 

5 

1.3 


1.26 X 10'° 

1.65 X 10 10 

6 

0.059 


1.10 X 10“ 

0.65 X 10 10 

7 

0.0011 

2.40 X 10“ 

0.26 X 10 10 


These relative densities are plotted in Figure 
3, together with the main nitrogen atmosphere 
and the Bates-Patterson hydrogen atmosphere 
[ Bates and Patterson, 1961]. The Bates-Pat- 
terson atmosphere, constructed for an exosphere 
temperature of 1500°K, is representative of the 
hydrogen distributions obtained by means of 
equation 1. The exponent k of the density scale 
is not specified because the absolute densities 
are not necessary here; k takes a different 
value for each curve. 


Discussion of Results 

The Monte Carlo curve represents diffusion 
through a slab medium bounded by vacuum on 
both sides. The density gradient in vacuum is 
zero (free particles), so that the density gradi- 
ent just within the medium must be zero also. 
This will hold for a distance into the medium 
that is of the order of a mean free path, i.e. 
until significant scattering occurs. This distance 
is only about 0.04 km at the lower boundary 
( 140 km) but about 25 km at the upper bound- 
ary (500 km) . The MC curve is based on seven 
widely separated points; therefore it does not 
show the change of slope that should occur 
within the border regions, but for accuracy it 
is terminated at 470 km. 

This density curve has a much steeper slope 
than that of Bates and Patterson. We first ask 
whether the net flow relation 

/ = DMnim/dz), ( 6 ) 

is satisfied for this curve, where Z) t is the diffusion 
coefficient at plane i, computed from the appli- 
cable formula of Hirschf elder et ah [1954, p. 539]; 


Z) 13 = 0.00268 


where 


[T 3 (M t + M 2 )/M, Af 2 ]‘ 


p<7u fija 


(l.l)*/; 


(7 


p = pressure in atmospheres. 

TV = kT/e lt . 

an, €u /k are potential parameters in angstroms 
and degrees Kelvin, respectively. 

is an integral given in Table 1 [Hirsch- 
feider et al, 1954, p. 1126] as a function of T. 

The product D i '(dn{H)/dz) i is given in Table 4. 
It is seen to decrease more and more rapidly with 
altitude, instead of remaining constant as we 
would normally expect. 

The explanation for this nonconstancy is 
believed to be the fact that the Chapman- 
Enskog kinetic theory, from which the diffusion 
coefficient is derived, becomes invalid when 
density and temperature changes are appreciable 
over a mean free path [Hirschj elder et al., 1954, 
pp. 18-20]. An effective diffusion coefficient in 
such a case does not appear to have been derived 
in the literature. The situation has been dealt 
with to some extent in the field of gas dynamics 
[Schaaf and Chambre, 1961], which we may 
qualitatively discuss to show its bearing on the 
atmospheric problem. 

The three main flow regimes of gas dynamics 
are the continuum flow regime, the free molecular 
flow regime, and the transition regime lying 
between them. The regimes are based on the 
range of values of a Knudsen number K, which 
is defined as the ratio of mean free path to some 
significant dimension of the flow field. Very 
roughly, the continuum regime corresponds to 
K <3C 1, the free molecule regime to K > 1, and 
the transition regime to intermediate values of 


TABLE 5 


X 

Mean Free Path 
l if km 

Scale Height 
hi, km 

* 

II 

j—'i 

1 

0.036 

25 

0.0014 

2 

0.10 

33 

0.0030 

3 

0.12 

38 

0.0031 

4 

0.21 

39 

0.0054 

5 

0.35 

41 

0.0085 

6 

2.0 

53 

0.038 

t 

25.5 

83 

0.31 


HYDROGEN DIFFUSION IN THE THERMOSPHERE 


163 


K. Though this division of gas dynamics pertains 
actually to flow past solid surfaces, it appears 
relevant in the present study to consider density 
scale height as a significant dimension. Knudsen 
numbers — l { /h are listed in Table 5, using 
mean free paths U taken from the computer 
calculation. K becomes appreciable compared 
with unity as altitude increases, so that condi- 
tions depart from those of the continuum regime 
and become those of the transition regime. The 
diffusion coefficient of the continuum regime 
will not then be valid. 

These observations do not constitute proof 
of the suggested explanation, of course, but only 
an argument for its plausibility. 

A direct measure of the departure of the dif- 
fusion coefficient formula from validity is af- 
forded by the entries of Table 4; since the 
product Di’dn(H) /dz decreases by a factor of 
11.5 between the altitudes of 140 and 296 km, 
D at 296 km is too small by this factor. 

As was noted earlier, gravity has been neg- 
lected. Its inclusion would decrease net flow, 
because molecular escape from the critical level 
near 500 km would then be an important factor 
in controlling the over-all escape rate. This 
factor is taken into account by Bates and Pat- 
terson, who assume a particular value for the 
molecularly limited escape flux. The slope of 
their curve is therefore correctly less than that 
of the MC curve at lower altitudes, where the 
effective diffusion coefficient is (approximately) 
the same for both. At higher altitudes, how- 
ever, inclusion of this factor would not decrease 
the slope of the MC curve to that of Bates and 
Patterson, because of the increasing difference 
in effective diffusion coefficient. The true curve 
should then have a slope intermediate between 
that of Bates and Patterson and of the present 
curve. 

The assumption adopted at the start, that 
the diffusing hydrogen is not in hydrostatic 
equilibrium, appears to play a lesser role in 
steepening the density curve than does the non- 
validity of the diffusion coefficient, though it is 
difficult to particularize its effect on the final 
result. 

An approximation was made in setting den- 
sity of the main atmosphere equal to zero above 
500 km. Collisions above this altitude would 
probably have a negligible effect on the diffusion 
process. 


Conclusions 

The Monte Carlo calculation of diffusion of 
hydrogen through a nitrogen atmosphere re- 
sults in a higher gradient of hydrogen density 
than current theories predict, even allowing for 
the fact that gravity was not included. The de- 
parture is ascribed to increasing non validity of 
the usual diffusion coefficient with altitude, be- 
cause the mean free path becomes significantly 
large compared with density scale height. The 
usual (calculated) diffusion coefficient is esti- 
mated to be an order of magnitude too low at 
300-km altitude. 

At low altitudes the true curve should have 
the slope of hydrogen distributions calculated 
in the usual ^ay, e.g. that of Bates and Patter- 
son. At higher altitudes, however, the true 
curve's slope will be intermediate between that 
of the present curve and that of Bates and 
Patterson. 

Since the diffusion-limited flow is found to 
be higher than in other calculations, the over- 
all escape flux of hydrogen will also be higher. 
This could account for the denser hydrogen 
envelope hypothesized by Donahue and Thomas 
[1963]. 
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